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A molecular theory of liquid water is identified and studied on the basis of computer simulation 
of the T1P3P model of liquid water. This theory would be exact for models of liquid water in 
which the intermolecular interactions vanish outside a finite spatial range, and therefore provides 
a precise analysis tool for investigating the effects of longer-ranged intermolecular interactions. We 
show how local order can be introduced through quasi-chemical theory. Long-ranged interactions 
are characterized generally by a conditional distribution of binding energies, and this formulation is 
interpreted as a regularization of the primitive statistical thermodynamic problem. These binding- 
energy distributions for liquid water are observed to be unimodal. The gaussian approximation 
proposed is remarkably successful in predicting the Gibbs free energy and the molar entropy of liquid 
water, as judged by comparison with numerically exact results. The remaining discrepancies are 
subtle quantitative problems that do have significant consequences for the thermodynamic properties 
that distinguish water from many other liquids. The basic subtlety of liquid water is found then in 
the competition of several effects which must be quantitatively balanced for realistic results. 



Introduction 

On the basis of the temperature dependence of its 
macroscopic properties, liquid water belongs to a class 
of liquids often referred to as "associated liquids." Fig. [I] 
gives one example of a characteristic distinction between 
liquid water and common organic solvents: the temper- 
ature dependence of the internal pressure. For van der 
Waals liquids, this temperature dependence is related di- 
rectly to the temperature dependence of the liquid den- 
sity. Organic solvents conform to this expectation, but 
liquid water exhibits a contrasting behavior, and that al- 
ternative behavior is attributed typically to the temper- 
ature dependence of water association through hydrogen 
bonding interactions. 

As an associated liquid, water is regarded broadly also 
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FIG. 1: The internal pressure of several solvents as a function 
of temperature along their vapor saturation curve.— For van 
der Waals liquids (du/dv) T ~ ap 2 . Organic solvents conform 
to this expectation, but water is qualitatively different. 



as a network liquid, though a definition of "network liq- 
uid" is usually not attempted. One aspect of the net- 
work view of liquid water is that models with defined 
short-ranged hydrogen bonding interactions that vanish 
precisely outside a finite range are designed specifically 
to study the network characteristics of that fluid. The 
work of Peery & Evansi gives a recent example that re- 
views and advances the forefront of statistical mechanical 
theory built explicitly on that network concept; this line 
of inquiry does indeed have an extended history, and the 
citations^^^ give further examples. The finite range 
of interactions is a common feature of lattice-gas models 
of liquid water 

The theory investigated here is motivated by a simple, 
surprising result which can be regarded as a theorem for 
models with intermolecular interactions that vanish out- 
side a radius A, whether or not the liquid under consid- 
eration would be regarded as a network liquid on some 
other basis. Specifically, it is that the excess chemical 
potential is precisely 

^ ox = /3 M c H x s (A)+lnp(n A =0) . (1) 

This excess chemical potential /i cx is the Gibbs free en- 
ergy per mole, in excess of the ideal contribution at the 
same density and temperature, of the one-component 
fluid considered. The thermodynamic temperature is 
T = l/k/3 where k is the Boltzmann's constant, /igg is the 
excess chemical potential for a hard-sphere solute, with 
A the distance of closest approach, at infinite dilution in 
the fluid. In the final term, p (n\ = 0) is the probability 
that a specific actual molecule in the fluid has zero (0) 
neighbors within the radius A. 

After articulation, the theorem Eq. [1] is straightfor- 
wardly understood, and a simple proof is given below. 
As an orienting comment, we note that Eq. [1] has been 
useful already in the theory of the hard-sphere fluids Of 
course, hard-sphere interactions are not involved literally 
in the interactions of physical systems, and it will tran- 
spire below that the hard-sphere contribution in Eq. [T] 
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simply serves to organize the statistical information, not 
as an assumed feature of an underlying intermolecular 
potential model. 

We emphasize that Eq. Q] is correct if intermolecular 
interactions vanish beyond the radius A, independently 
of whether the liquid under consideration is a network 
liquid. Real network liquids have strong, short-ranged 
binding interactions, but invariably longer-ranged inter- 
actions as well. This paper investigates the simplest ap- 
proximate extension of the Eq. [T] to real liquids, an ex- 
tension which is realizable, and which is conceived on the 
basis of observable data, not from consideration of an as- 
sumed interaction model. The specific issues we address 
quantitatively include whether the effects of longer-than- 
near-neighbor interactions are simple for realistic models 
of liquid water, and whether reasonable values of A can 
be found that make the statistical theory of liquid water 
particularly simple with realistic intermolecular interac- 
tions. 



Statistical Thermodynamic Theory 

We seek /i° x on the basis of accessible simulation data. 
Since Eq. [1] contains we consider evaluating that 

quantity on the basis of simulation of the realistically 
modeled liquid of interest. For that purpose, the poten- 
tial distribution theorem (PDT ) 10 ' 11 difference formula 
yields 
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M ° X) = P (n x = 0) / P (e\n x = 0) e^de . (2) 



If we take /xgg as known,— Eq. ([2]) is an equation for 
/i cx . P (e\ri\ = 0) is the conditional probability of the 
binding energy of a specific molecule for the event that 
nx = 0, i.e., the specific molecule considered has zero 
(0) neighbors within radius A. p (n\ — 0) is the marginal 
probability of that event. If the intermolecular interac- 
tions vanish for ranges as large as A, those binding en- 
ergies are always zero (0). This yields Eq. Q] under the 
assumptions noted. 

For general interactions, we regard the conditioning 
of Eq. [5] when A > as a regularization of the statistical 
problem embodied in Eq. ([2]) when A = 0, which is impos- 
sible on the basis of a direct calculation. After regulariza- 
tion, the statistical problem becomes merely difficult. For 
A — > oo, a gaussian distribution model for P (e\n\ = 0) 
should be accurate since then many solution elements will 
make small, weakly-correlated contributions to e. The 
marginal probability p (n\ = 0) becomes increasingly dif- 
ficult to evaluate as A becomes large, however. For A 
on the order of molecular length scales typical of dense 
liquids, a simple gaussian model would accept some ap- 
proximation error as the price for manageable statistical 
error. If P (e\n\ = 0) is modeled by a gaussian of mean 
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We regard the formulation Eq. [3] merely as a parsimo- 
nious use of statistical information that might be ob- 
tained from a simulation record. This simple model mo- 
tivates the following analyses. 

An essentially thermodynamic derivation of Eq. [21 one 
that illustrates the connection with quasi-chemical the- 
ory, was given previously. 13 Our result there was formu- 
lated as 
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where a indicates the distinguished molecule and 'w' 
the molecules of the solvent medium, for example water. 
Comparison of Eqs. [5] and 2] then suggests the identifica- 
tion 
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Thus the marginal probability of the event (n\ = 0) di- 
rectly interrogates chemical contributions involving the 
chemical equilibrium ratios K^s on the basis of the as- 
sumed forcefield. This correspondence is indeed a ba- 
sic result of constructive derivations of quasi-chemical 
theory— Given adequate simulation data, explicit evalu- 
ation of individual K^s is not required. But more basic 
chemical analysis of the if n 's, e.g., using electronic struc- 
ture method ^l^j^^w^w^Oj^i^^s j g not implemented ei- 
ther. We add for completeness that previous discussions 
denoted p(n\ = 0) = Xq, with x n oc K n p w n the proba- 
bility that the distinguished molecule has a coordination 
number of n, and K =l*MMM££L 

These equations correctly suggest that this formulation 
is fully general for the classical statistical thermodynamic 
problem considered. This formulation doesn't make a lat- 
tice structural idealization, doesn't assume that the inter- 
action potential energy is pairwise decomposable, doesn't 
assume that the interaction potential energy model pro- 
vides a handy definition of 'H-bonded,' and the coordina- 
tion definition — the occupancy within a radius A — is 
observational on the basis of molecular structure. Effec- 
tive values of A should be established eventually on the 
basis of the statistical information secured. 

Of course, the chemical potential provides the partial 
molar entropy 
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On the right side the quantities (e) , and p (the fluid pres- 
sure) are mechanical properties which are directly avail- 
able from a simulation record. 

It is an important but a tangential point that the 
distributions P (e\n\ = 0) to which this theory directs 
our attention are different from the p(v) that have been 
frequently considered following the work of Rahman & 
Stillinger^ see also the review of StillingerJS p(v) is the 
distribution with the interpretation that Jjj 2 p(v)dv = 
n(v2,vi) is the number of molecules neighboring a par- 
ticular one with pair interaction v in the range (v2,i>i). 
p(v) has been used to convey reasonable definitions of 
'H-bonded' on the basis of simulation models of water, 
and it is well-designed for the purpose. It doesn't supply 
the correlations that the distributions P (e\n\ = 0) use 
to determine the free energy. 



Simulation Data 

To test these ideas, 50000 snapshots from 
simulations 2 -^ performed at 298, 350 and 400 K 
were analyzed. Calculations of and p(n\ = 0) were 
carried out as described by Paliwal, et alM- The inter- 
action energy of a water molecule was calculated as the 
sum of pairwise additive van der Waals and electrostatic 
energies. Lennard-Jones interactions were considered 
within 13 A, and standard long range corrections were 
applied beyond this distance. Electrostatic interactions 
were evaluated using the Ewald summation method. 

Statistical uncertainties in /i ox were computed by di- 
viding the total number of snapshots into five blocks and 
evaluating block averages. For comparison, /^ ox was de- 
termined from the histogram overlap method after eval- 
uating also the uncoupled binding-energy distribution 
p(°) which was obtained by placing 5000 randomly 

oriented water molecules at randomly chosen locations 
in each snapshot. Table Q] collects the individual terms 
for the gaussian model, Eq. ([3]), at each temperature, 
and additionally the entropy evaluated from Eq. [7J The 
observed dependence on A of the free energy at each tem- 
perature is shown in Fig [3] 



Results 

Fig. [5] shows that the unconditioned distribution P (e) 
displays positive skew, but the conditioning diminishes 
that skew perceptibly, as expected. P {e\n\ = 0) is least 
skewed for the largest A. All of these results are quali- 
tatively unlike the binding-energy distributions obtained 
for a two dimensional model of liquid water— which has 
been studied furtheri 30 i 31 That previous result is multi- 
modal, behavior not seen here. 

The distributions P (e\n\ — 0) shown in Fig. [5] are 
also qualitatively unlike the pair interaction distributions 
j?(i>) 25 ' 26 which are different, as was discussed above. 
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FIG. 2: Probability density P(e\n\ = 0) of the binding en- 
ergy of a water molecule in liquid water at 298 K. Se = 
e — (e\ri\ = 0), A = 0.0, . . . 0.30 nm, from top to bottom with 
successive results shifted incrementally downward by 2 for 
clarity. The solid lines are the gaussian model for each data 
set. 



The conditioning affects both the high-e and low-£ 
tails of these distributions. The mean binding energy 
(e\n\ = 0) increases with increasing A [Tabled], so we 
conclude that the conditioning eliminates low-e configu- 
rations more than high-e configurations that reflect less 
favorable interactions. Nevertheless, because of the ex- 
ponential weighting of the integrand of Eq. ((2]) and the 
large variances, the high-e side of the distributions is 
overwhelmingly the more significant in this free energy 
prediction. 

The fluctuation contribution exhibits a broad maxi- 
mum for A < 0.29 nm, after which this contribution de- 
creases steadily with increasing A [Table |J . Evidently 
water molecules closest to the distinguished molecule, 
i.e., those closer than the principal maximum of oxygen- 
oxygen radial distribution function, don't contribute im- 
portantly to the fluctuations. This is consistent with 
a quasi-chemical picture in which a water molecule and 
its nearest neighbors have a definite structural integrity. 
Fig. 2] shows x n oc K n p w n ] the most probable coordina- 
tion number is either 2 or 3 when A=0.30 nm, but it is 4 
when A=0.33 nm. Nevertheless, the breadth of this dis- 
tribution is remarkable. The coordination number n=6 
is more probable than n=2 for the A = 0.33 nm results 
of Fig. H 

The magnitude of the individual contributions to /i cx 
are of the same order as the net free energy; the mean 
binding energies are larger than that, as are the variance 
contributions in some cases. The variance contributions 
are about half as large as the mean binding energies, with 
opposite sign. It is remarkable and significant, therefore, 
that the net free energies at all temperatures are within 
roughly 12% of the numerically exact value computed 
by the histogram-overlap method. The predicted excess 
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TABLE I: Free energy contributions in kcal/mol associated with the gaussian model. The bottom value of each temperature 
set gives the corresponding free energy evaluated by the histogram overlap method. The rightmost column is the excess entropy 
per particle. Note that the entropy here raises the net free energy (see Eq. [7)l , and the magnitude of the entropy contribution 
TS e */n ranges here from 60-75% of the magnitude of the free energy /i cx . 
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0.2600 


2.93 


-0.04 


-19.74 


9.87 


-6.98 ± 0.02 


-5.88 




0.2650 


3.15 


-0.13 


-19.68 


9.93 


-6.73 ± 0.03 


-6.30 




0.2675 


3.25 


-0.20 


-19.59 


9.97 


-6.57 ± 0.04 


-6.57 




0.2700 


3.35 


-0.31 


-19.46 


9.98 


-6.44 ± 0.03 


-6.79 




0.2725 


3.45 


-0.44 


-19.27 


9.97 


-6.29 ± 0.03 


-7.05 




0.2750 


3.56 


-0.60 


-19.03 


9.92 


-6.15 ± 0.03 


-7.28 




0.2775 


3.67 


-0.78 


-18.73 


9.83 


-6.01 ± 0.04 


-7.52 




0.2800 


3.78 


-0.98 


-18.39 


9.71 


-5.88 ± 0.03 


-7.74 




0.2900 


4.23 


-1.96 


-16.75 


8.89 


-5.59 ± 0.04 


-8.23 




0.3000 


4.71 


-3.09 


-14.93 


7.77 


-5.54 ± 0.06 


-8.31 




0.3100 


5.22 


-4.27 


-13.21 


6.67 


-5.59 ± 0.18 


-8.23 




0.3200 


5.76 


-5.45 


-11.65 


5.54 


-5.80 ± 0.37 


-7.87 




0.3300 


6.32 


-6.64 


-10.30 


4.78 


-5.84 ± 0.97 
-6.49 


-7.81 
-6.71 


350 


0.2600 


3.06 


-0.05 


-18.43 


9.44 


-5.98 ± 0.02 


-5.65 




0.2625 


3.15 


-0.09 


-18.41 


9.46 


-5.89 ± 0.02 


-5.78 




0.2700 


3.45 


-0.33 


-18.14 


9.47 


-5.55 ± 0.02 


-6.27 




0.2750 


3.66 


-0.62 


-17.73 


9.35 


-5.34 ± 0.02 


-6.57 




0.2800 


3.87 


-0.99 


-17.15 


9.10 


-5.17 ± 0.01 


-6.82 




0.2900 


4.32 


-1.92 


-15.67 


8.24 


-5.03 ± 0.04 


-7.02 




0.3000 


4.80 


-3.00 


-14.05 


7.20 


-5.05 ± 0.06 


-6.99 




0.3100 


5.29 


-4.13 


-12.48 


6.14 


-5.18 ± 0.10 


-6.80 




0.3200 


5.81 


-5.28 


-11.01 


5.24 


-5.24 ± 0.22 


-6.72 




0.3300 


6.36 


-6.41 


-9.74 


4.47 


-5.32 ± 0.45 
-5.83 


-6.60 
-5.87 


400 


0.2600 


3.03 


-0.06 


-17.19 


8.97 


-5.25 ± 0.02 


-5.21 




0.2625 


3.13 


-0.10 


-17.16 


8.99 


-5.14 ± 0.03 


-5.35 




0.2700 


3.41 


-0.35 


-16.89 


8.95 


-4.88 ± 0.02 


-5.67 




0.2750 


3.61 


-0.63 


-16.49 


8.80 


-4.71 ± 0.03 


-5.89 




0.2800 


3.81 


-0.98 


-15.96 


8.52 


-4.61 ± 0.03 


-6.01 




0.2900 


4.23 


-1.87 


-14.60 


7.69 


-4.55 ± 0.05 


-6.09 




0.3000 


4.68 


-2.89 


-13.10 


6.72 


-4.59 ± 0.05 


-6.04 




0.3100 


5.14 


-3.97 


-11.64 


5.78 


-4.69 ± 0.06 


-5.91 




0.3200 


5.62 


-5.06 


-10.31 


4.95 


-4.80 ± 0.14 


-5.77 




0.3300 


6.13 


-6.12 


-9.12 


4.27 


-4.84 ± 0.18 
-5.31 


-5.72 
-5.13 



entropy, S CVL /nk, (Fig. [3]) is in error by 17% at the low- 
est temperature considered and more accurate than that 
at the higher temperatures. The constant-pressure heat 
capacity predicted by the model is roughly 30% larger 
than the result obtained using the numerically exact val- 
ues for S ex /nk in the same finite-difference estimate. A 
mean-field-like approximation that neglects fluctuations 
produces useless results. 

The consistent combination of the various terms is an 
important observation. For example, the packing con- 
tribution, /igg (A), is sensitive to the interaction poten- 
tial energy model used, and primarily through differ- 
ences in the densities of different models at the higher 
temperatures.— We found that use of values of /i£f s (A) 
corresponding to a more accurate model of liquid water 12 
spoils its consistency with the other contributions evalu- 
ated here with the TIP3P model, and that appreciably 



degrades the accuracy of the whole. 

If kT\np(n\ — 0) is zero [TableQ], the hard-sphere con- 
tribution //gg (A) is ill-defined. As a general matter, the 
sum /ugg (A) + kT\up(n\ = 0) cannot be identified as a 
hard-sphere contribution. Since these terms have oppo- 
site signs, the net value can be zero or negative, and 
those possibilities are realized in Tabic fl] To define the 
hard-sphere contribution more generally, we require that 
Mhs M ^ e continuous as p(n\ = 0) — > 1 with decreasing 
A. All other terms of Eq. ([3]) will be independent of A 
for values smaller than that, and it is natural to require 
that of j«gg (A) also. 

From Fig. [3J we see that A > 0.30 nm identifies a larger- 
size regime where the variation of the free energy with 
A is not statistically significant. Although we anticipate 
a decay toward the numerically exact value for A — > oo, 
the statistical errors become unmanageable for values of 
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A much larger than 0.30 nm. When A = 0.30 nm a sig- 
nificant skew in P(e\n\ = 0) is not observed, as already 
noted with Fig. [2] 

The predicted free energy fj, cx is then distinctly above 
the numerically exact value, suggesting that the gaussian 
model predicts too much weight in the high-e tail. We 
further examined the high-e tail of P (e\n\ = 0) for larger 
values of A by carrying-out Monte Carlo calculations un- 
der the same conditions but with one water molecule car- 
rying an O-centered sphere that prohibits other O-atoms 
within a radius A. Those results are shown in Fig. [5] The 
accuracy of a gaussian model P (e\n x = 0) is clearly born- 
out, but the gaussian model predicts slightly too much 
probability in the extreme high-e region particularly for 
A = 0.33 nm. The radial distribution of outer-shell neigh- 
bors in the conditioned sample (Fig. [6]) shows that in that 
case the mean water-molecule-density is structured as a 
thin radial shell with roughly n ~ 15 near-neighbors. 

The design strategy here is that the theory should be 
correct in the A — > oo limit. The accuracy of this ap- 
proximate theory should then be judged for the largest 
values of A that are practical. Consulting the direct in- 
vestigation of the high-e tail (Fig. [5]), we see that the 
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FIG. 4: Coordination number probabilities, x n = p{n\ = n). 
Note that for A=0.30 nm, the most probable coordination 
number is 2 or 3.—^ For A=0.33 nm, the most probable 
coordination number is n=4. But n=4 and n~5 are nearly 
equally probable at T—298 K, and n = 6 is then 17% of the 
whole. With the TIP3P model, as T is lowered from 400 K 
to 298 K there is a clear tendency to enhance probability of 
higher coordination number cases. 
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FIG. 3: upper panel: Dependence of the free energy ^t ex pre- 
dicted by the gaussian model on the conditioning radius A. 
The horizontal lines are the numerically exact results. The 
shaded areas indicate approximate 95% confidenc 
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probability density drops off by roughly e -10 from the 
mode to the highest-e data there. But the integrand fac- 
tor of e° s of Eq. [2] grows by roughly e +20 across that 
range; these distributions are remarkably wide even after 
the conditioning. Thus a statistical model that extends 
those data is still necessary. 

In the context of implementations of quasi-chemical 
theory, it has been argued that the cut-off parameter A 
may be determined on a variational basis. Inner- and 
outer-shell contributions have countervailing effects, so 
where the net-results are insensitive to local changes 
of A distinct approximations may be considered well 
balanced^ Here, the contributions corresponding to the 
inner-shell contribution in quasi-chemical theory are ob- 
tained numerically exactly so that argument isn't com- 
pelling. 

A previous quasi-chemical analysis of liquid water^ 
which utilized ab initio molecular dynamics (AIMD) 
data, provides an interesting comparison with the present 
results. In initiating analysis of AIMD, that previous ef- 
fort addressed a much more complicated case. Several 
contributions had to be estimated on the basis of slightly 
extraneous input, e.g., dielectric model calculations and 
rough estimates of the effects of outer-shell dispersion 
contributions. Nevertheless, there was some commonal- 
ity in that both mean-field and fluctuation (dielectric) ef- 
fects were involved in the outer-shell contributions. The 
net accuracy in that previous case was similar to that 
observed for the simpler case, more thoroughly analyzed 
here. The evaluated free energy was too high but by less 
than 1 kcal/mol at moderate temperatures. 

Review of the theory Eq. [3] identifies a further remark- 
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FIG. 5: The fine black points are estimates of P (e\n\ — 0) 
obtained from Monte Carlo simulation of TIP3P water with 
one water molecule centering a hard-sphere which perfectly 
repels other water oxygens from a sphere of radius A. (up- 
per panel) A = 0.30 nm; (lower panel) A = 0.33 nm. The 
larger composite dots were obtained directly from molecular 
dynamics simulations with no constraint. The solid curve is 
the gaussian model for the expanded Monte Carlo data set, 
and confirms both the overall accuracy of the gaussian model 
and the estimated conditional mean and variance of Table I. 
Note that the data are slightly less than the gaussian model 
in the high-£ tail. T=298K. 



able observation. If we follow that logic in reverse, ex- 
ploiting the fact that /i cx is an experimentally known 
quantity, then Eq. [3] provides a physically transparent 
theory for /igg (A) the object of central interest in theo- 
ries of hydrophobic effects i 12 ' 34 As a theory of hydropho- 
bic effects, Eq.[3]offers a specific accounting for the effects 
of the distinctive coordination characteristics of a water 
molecule in liquid water, and a specific accounting of the 
role of longer-than-near-neighbor-ranged interactions. 




FIG. 6: Distribution of water oxygen atoms radially from a 
distinguished water oxygen atom; evaluated with the sample 
corresponding to the condition n\ — 0, with A = 0.30 nm, 
for O-atom neighbors of the distinguished O-atom. The 
black, dotted curve is the running coordination number for 
the T=298K case. 



Conclusions 

Viewed from the perspective of Eq.[TJ which is precisely 
correct for liquid water models with interactions shorter- 
ranged that A, the quantitative contribution from longer- 
ranged interactions for the TIP3P model is essentially 
identical to the experimental value of /x ex for reasonable 
values of A ~ 0.33 nm, independently of temperature. 
If these longer-ranged contributions were modeled as a 
mean- field contribution {e\n\ — 0) only, the error would 
be nearly 100%. On the basis of the present work, it is 
evident that Ref. 35, in adopting the value A = 0.28 nm 
for /i^ g (A), assumed that quantitative chemical contribu- 
tions reflecting short-ranged local order were negligible, 
independently of the local environment. This work then 
clarifies that for bulk water A — 0.26 nm, rather than 
0.28 nm, if those chemical contributions are to be negli- 
gible. Moreover, suitable values of A require additional 
scrutiny for different local environments. 

The present theory helps to gauge the significance of 
the various features of the interactions for the molecular 
of liquid water. As an example consider the case of an 
interaction model with near-neighbor-ranged interactions 
only. Then Eq.[T] applies. Suppose that the parameters of 
the model are adjusted to match precisely the observed 
coordination number probabilities (Fig. [4]) for realistic 
models. If the packing contribution is obtained from 
physical dat a 32 i 35 or from simulation of realistic models 
of liquid water^ then it will have positive magnitudes 
consistent with the values of Table |TJ Then we must an- 
ticipate that the near-neighbor-ranged interaction model 
will yield errors in the free energy of roughly 100%. Con- 
trapositively, if the parameters of such a model are ad- 



7 



justed to yield an accurate free energies, /i cx and /xgg (A), 
the results above together with the theorem Eq.Q] suggest 
that the predicted coordination number distribution - 
specifically p (n\ = 0) = xq — should differ significantly 
from results for models with realistic longer-ranged inter- 
actions. The implication of the analysis presented here 
for network models of liquid water is that a detailed treat- 
ment of local order alone - e.g., invoking specific short- 
ranged hydrogen-bonding interactions that capture geo- 
metric features of the tetrahedral structure of ice but not 
involving longer-ranged interactions - does not generally 
describe the range of thermodynamic properties that dis- 
tinguishes water from organic solvents. 

A more specific conclusion is that the observed 
P (e\n\ = 0) are unimodal. The gaussian model inves- 
tigated here is surprisingly realistic; the molar excess 
entropies are predicted with an accuracy of magnitude 
roughly k (Boltzmann's constant). The predicted free 
energies are more accurate than 1 kcal/mol. The distinc- 



tions from gaussian form are subtle quantitative prob- 
lems that have significant consequences for the thermody- 
namic properties that distinguish water from many other 
liquids. Those quantitative subtleties remain to be fur- 
ther analyzed. 
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